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ABSTRACT 

The structure of accretion discs around magnetic T Tauri stars is numer¬ 
ically calculated using a particle hydrodynamical code, in which magnetic 
interaction is included in the framework of King’s dimagnetic blob accretion 
model. Setting up the calculation so as to simulate the density structure of a 
quasi-steady disc in the equatorial plane of a T Tauri star, we find that the 
central star’s magnetic field typically produces a central hole in the disc and 
spreads out the surface density distribution. We argue that this result suggets 
a promising mechanism for explaining the unusual flatness (IR excess) of T 
Tauri accretion disc spectra. 

Key words: accretion, accretion discs - stars: pre-main-sequence - stars: 
magnetic fields. 


1 INTRODUCTION 

It is now quite widely accepted that in the formation process of low-mass (M < 2Mq) stars, 
a disc accretion phase is typically present (see Hartmann 1998 for an extensive review and 
references). In particular, the classical T Tauri stars (CTTS) are thought to be still sur¬ 
rounded by discs and the emission properties of some typical objects suggest that these 
discs are not merely “passive” dust bodies, irradiated by the central star, but are rather 
bona-fide energy producing, viscous accretion discs, similar to the ones present in mass 
transferring close binary systems, and whose modeling dates back to the classical works of 
Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974). Because of their ubiquity and 
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of some interesting theoretical challenges they pose, accretion discs have remained at the 
centre stage of astrophysical research (see e.g. Frank, King & Raine 1992; Papaloizou & Lin 
1995; Lin & Papaloizou 1996 for reviews). 

One particular issue in this context is the fact that in many CTTS the measured spectrum 
presents an unusually high emission in the infra red (called IR excess), far more than expected 
from disc reprocessing of stellar light or viscous heating of the disc (see Bertout 1989, 
Beckwith et al. 1990, Hartmann et al. 1994). Our wish to investigate this problem, as well as 
some other facets of CTTS discs, has been the main motivation for the work reported on in 
this paper. It consists of setting up a numerical calculation of the properties of an accretion 
disc around magnetic CTTS, using a scheme of magnetic interaction, originally proposed in 
the context of diamagnetic accretion by King (|1993|) . 

This description of the magnetic interaction assumes that as material moves through the 
magnetosphere it interacts with the local magnetic held via a velocity dependent acceleration 
(force per unit mass, /mag) of the form: 


/mag = -K[U- U/]_L, 


( 1 ) 


where u and Uf are the velocities of the material and magnetic held lines respectively, iC is a 
suitable ’’magnetic drag” coefficient (see below) and the suffix T refers to the vector compo¬ 
nent perpendicular to the held line. The magnetic acceleration, as expressed in equation |^, 
is intended to represent the dominant term of the magnetic interaction, with K containing 
the relevant parameters determining the ehective magnetic time-scale. 

Within the above description (HD there are still a number of diherent possibilities to model 
the inner disc - magnetosphere interaction. In this paper we shall use the diamagnetic blob 
accretion (DBA) model, although it is possible to formulate the diametrically opposite case 
(complete magnetic penetration of the disc) in the general form (|^) as well. Wynn, Leach & 
King (pOOlD show how to calculate the appropriate coefficient K for the latter case and we 
plan to repeat our calculations in the future, using this prescription. 

In the DBA approach one assumes that the fluid constituting the accretion flow, i.e. 
the disc and its surroundings, is composed of blobs immersed in a dilute interblob plasma 
(see Aly & Kuijpers 1990). The blobs behave diamagnetically in the presence of the stellar 
magnetic held and thus suffer a surface drag force acting on them. The model has since 
its introduction been applied to various systems, notably to the intermediate polar class 
of cataclysmic variables (CV). Wynn & King (|1995|) and Wynn, King & Horne ( 1997 ) 
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incorporated the diamagnetic drag force into a particle hydrodynamical numerical code 
(HYDISC), originally developed by Whitehurst ( |1988|) to simulate the accretion disc in non¬ 
magnetic close binary systems. They found important properties of the white dwarf’s spin 
evolution and its effect on the disc and explicitly applied their hndings to the moderately 
magnetic CV AE Aqr. A more recent study using this approach was its application to the 
system EX Hya ([King fc: Wynn 19991) . 

The DBA approach was extended to the study of accreting T Tauri stars for the hrst 
time by King & Regev (|1994| ), hereafter KR. They performed calculations of individual h\oh 
(that is, ballistic) orbits, including the interaction with magnetic loops of the central star 
and have shown that this mechanism can eject blobs from the system, in directions pointing 
away from the disc plane. By estimating the angular momentum thus expelled and including 
it in the overall angular momentum budget KR found that a stellar spin equilibrium value 
is compatible with the observations (i.e. ~ an order of magnitude less than the breakup 
value) if the magnetic loops extend to a few stellar radii. Pearson & King (|1995| ), hereafter 
PK, generalised the above work into a A^-particle simulation, in a similar manner to the one 
mentioned above for CVs. This work also focused on the issue of the slow observed T Tauri 
spin rates and conhrmed the idea that such equilibrium spin rates can be achieved by the 
expulsion of matter from the disc until the corotation radius (see below) coincides with the 
edge of the magnetic loop. The additional new Ending, not anticipated by KR, of this work 
was that the ejection of material comes to a halt as the star approaches its spin equilibrium 
value. 

The problem of the slow spin rates and rotational evolution of T Tauri stars as a result 
of the star’s magnetic interaction with its accretion disc has been also approached, in a 
number of works, by analytical and semi-analytical methods, that is, stopping short of 
multi-dimensional numerical simulations. Konigl (1991), Cameron & Campbell (1993) and 
Armitage & Clarke (1996) used quite different approaches and all found that a stellar dipole 
magnetic held of strength < 1 kC is able to regulate the stellar spin to a quasi-static value. 


in the observed range. In the process of the magnetic interaction the inner accretion disc 
gets disrupted and the inner (abrupt) termination radius of the disc was estimated in the 
above papers. 

KR and PK used the DBA approach focusing on the ejection of mass from the disc, 
in an effort to link the low spin rates of T Tauri stars with outhows from young stellar 
objects (YSO) within a unihed scenario. In the present work we utilise the DBA model in 
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the calculation of the properties of the steady (or quasi-steady) accretion disc itself, that 
is, the material that remains close to the equatorial plane while spiraling in due to viscous 
torques, after t\ie spin period has already stabilised. Some pecularities of T Tauri light-curves 
(i.e. the temporal luminosity variations) have already been treated using a model based on 
single blob ballistic orbits, which leave the the disc plane but ultimately return and impinge 
on the disc surface ([Ultchin et ah 1997|) . In order to investigate whether the DBA model 
is consistent with some of the basic properties of T Tauri spectra (see Bertout 1989 and 
references therein) and in particular the above mentioned IR excess, we have performed 
numerical simulations using a code based on HYDISC, suitably modihed so as to adapt it 
for the problem at hand. 

This paper is organised as follows. In the next section we discuss the basics of blob 
dynamics and estimate the relevant time-scales. §3 describes the numerical code used in the 
accretion flow simulation and the procedure for hnding the spectrum emitted by such flows. 
In §4 the results of our simulations for a score of parameter values are described in some 
detail and hnally, §5 summarises this work in comparison to the results of other approaches. 


2 DIAMAGNETIC ACCRETION FLOW DYNAMICS 

Following KR and PK we model the material of a T Tauri disc as being composed of gas blobs 
immersed in a tenuous interblob plasma. These blobs are assumed to behave diamagnetically 
in the presence of the magnetic held of the central star, which threads the plasma. Drell, 
Folley & Ruderman (1965) showed that the typical time-scale on which a diamagnetic object 
of mass m (i.e. a blob) loses energy, when it moves in a magnetic held B, is 

CAfn 

where ca = \J/ (4vrpp) is the Alfven speed of the interblob plasma, whose density is Pp, 
and lb is the blob size. There is a minimum condition on the conductivity for this equation 
to hold and, as it was shown by KR, this condition is satished for all reasonable parameters 
in the conditions appropriate for T Tauri accretion. 

This allows to cast the diamagnetic force per unit mass, acting on a blob, into one simple 
formula, whose form is in accordance with equation |^: 

fd = -Kd [Ufe - U/]_L h. 


( 3 ) 
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where = 1/td (see below) is the diamagnetic drag coefficient. The magnetic drag force 
obviously vanishes if the blob velocity u;, relative to the velocity of the local held line uj is 
parallel to the held lines, it is proportional to the perpendicular (to the local held) component 
of this relative velocity and is directed in the normal direction h to the held. This is the 
meaning of the shorthand vector notation in equation ^ Its explicit component form will be 
given below. 

Assuming that the magnetic held rotates with the star and does not change on the 
time-scale of interest, the above equation is best explicitly written in the frame of reference 
rotating with the central star and had the form 

frf =-A'd [v - (v ■ b)b], (4) 

where v is the blob velocity in the rotating frame, that is, relatively to the held lines, and b 
is a unit vector in the local held direction (see KR). 

Substituting explicitly the above form of ca and m = llpb, where pb is the blob density, 
in the drag coefficient one gets 

Rd(x) = 2^py^pblb)-^ B = Kffs), (5) 

where the second equality parametrizes in terms of a constant coefficient K and a spatially 
dependent (x denotes the position vector) function, which is of order 1. The blob parameters 
and their space dependence, as well as the plasma density, are rather uncertain. The only 
variable which can be explicitly written is the magnetic held (in the approximation that it is 
hxed in the rotating frame, that is, anchored to the star), because it can be chosen a priori. 

In this work we shall mainly use a dipolar magnetic held conhguration and group, for 
the sake of convenience, all the other variables into a single hxed parameter, incorporating 
in it most of our ignorance of the complex plasma instabilities, which are supposed to create 
and shape the blobs. Within this framework we have 

Kd = K{R/R,)-^[1 + ^{z/R)Y/\ ( 6 ) 

in cylindrical coordinates {R, z) whose z axis coincides with the dipole axis. R* is the stellar 
radius and 

K = V^fi'P(nh)-'Ba, (7) 

where Bo is the polar held strength. This K is assumed to be constant and will serve as 
a parameter of our problem. For a discussion of the diherent options to chose see e.g. 
Wynn, King & Horne (1997). 
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We focus in this work on the case of an aligned dipole, that is, when the dipole axis 
coincides with the disc axis, but will briefly discuss (in §4.3) the inclined dipole case as well. 

The magnetic drag term has to be added to the blob equations of motion, which in¬ 
clude also gravitational, centrifugal and coriolis terms (in the rotating frame). Scaling the 
independent variables by their natural values, that is, lengths by the corotation radius, 
Rco = (with M* the central star’s mass and Vl its rotational angular velocity) 

and time by 1/ff, the following blob equations of motion are found ([Ultchin et al. 19^) : 


X 

x = 2y + x-^- fcrf(x) [i: - (x ■ b) b^, 
y = -2x + y- ^ - ^d(x) [y- (x- b) by], 
Z = - fcd(x) [i - (x ■ b) b^], 


( 8 ) 

( 9 ) 

( 10 ) 


These constitute the nondimensional equations of motion for the cartesian components, 
in the rotating frame, of the blob trajectory x(f) and are thus ready to be incorporated into 
HYDISC (see in the next section). R, the cylindrical radius, is obviously R = {x‘^ -|- 
and the calculation of b is quite simple for an aligned dipole. For a tilted dipole, a somewhat 
more involved (but straight forward) calculation is needed. 

The drag term has been written as kd, the nondimensional version of Kd (of equation 
The explicit form of this function, following from equations ^ and |^, brings out the physical 
meaning of the constant coefficient k appearing in front of the nondimensional spatially 
dependent part of the magnetic held 

kd{x) = k {R/RY-^ [1 + 3(^/i^)]'/^ (11) 

which is 


k = 2.82 




Pp 


1/2 


hpb 


-1 


lOOG/ \l0“^°gcm“3 


100 g cm 


-2 


^rot 


106 s; Wg(i?*,iooG)' 


( 12 ) 


Here reasonable (order of magnitude) estimates of the parameters for accreting T Tauri 
systems have been used to scale the relevant variables. P* ~ pot is the stellar spin period and 
Bq the polar magnetic held strength. The disc density scaling value is estimated from models 
of standard accretion discs around T Tauri stars, using a representative accretion rates of 
M ~ 5 lO“®M 0 yr“^ ( jBasri fc Bertout 1989|) . The blob parameters are rather uncertain, but 
it is reasonable that the blob size is probably of the order of the vertical scale height in the 
disc while its density should be several orders of magnitude larger than that of the disc. 
This gives the value with which we scale lb Pb- In any case, our results will be seen to depend 
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only on the parameter k (and several other quantities related to the numerical model) and 


hence on the particular combination given in equation 

Time-scale estimates can now be easily seen to demonstrate that for fc ~ 1 the magnetic 
drag time-scale near the star (where i? ~ i?* and thus the contribution of the space depen¬ 
dent part, arising from eq. 0, is very close to 1) for a polar held of ~ lOOG, denoted here 
by tmag(-R*, lOOG), is much larger than the dynamical time-scale there, since the latter is 
much shorter than Got (T Tauri stars are slow rotators). Thus we do not expect a signihcant 
effect near the star due to the magnetic drag in this case. Since the magnetic held decays 
(and thus the magnetic drag time scale grows) faster with R than the dynamical time-scale, 
the ehect of the magnetic held on the motion will be small everywhere for fc ~ 1 

Let us examine now what happens at the corotation radius. From the functional depen¬ 
dence of the dipole held (eq. 0) we have the term [Rco/Rfij ^ = (G/G^)^, with the ratio of 
the G-s evaluated at the stellar surface. This term is of the order of 10“^ in T Tauri stars, 
but trot; by dehnition, is equal to the dynamical time scale at the corotation radius. If one 
chooses k ~ 100, a signihcant ehect of the magnetic force at corotation is expected and 
obviously more so inside corotation. Thus we expect a signihcant modihcation of the inner 
accretion disc in this case. 

These time-scale estimates provide the rationale for our choice of the k parameter range 
(or Bq with all the other parameters constant) in the simulations (described in §4). Note 
that the viscous time-scale did not enter our considerations here. In accretion disc theory 
it is always assumed that this time-scale is much larger than the dynamical time-scale and 
consequently if the the magnetic time-scale approaches the dynamical one they are both 
much shorter than the viscous time scale. However, in regions where the magnetic time- 
scale is extremely long (far out in the disc, or even quite close in, if the held is weak enough) 
it might be even longer than the viscous time-scale and thus the magnetic drag would not 
interfere with the viscous evolution of the disc. 


3 THE NUMERICAL METHOD 
3.1 General considerations 


To facilitate an accretion how simulation in our system we adapt HYDISG, a numerical par¬ 
ticle code originally written by Whitehurst (|198^) to study accretion in GV, to the problem 
of a single central magnetic accretor (a T Tauri star). The introduction of the magnetic 
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interaction into the particles’ equations of motion, using the DBA approach, has already 
been explained in the previous section. Our blobs are being considered as the particles of 
the particle scheme, utilized in HYDISC to simulate the hydrodynamical aspects of the flow, 
like pressure and viscosity. As it was already mentioned in the Introduction, similar schemes 
based on HYDISC have been designed for simulating accretion flows in magnetic CV (|W ynn 


fc King 19951; [Wynn, King fc Horne 1997|) and also for investigating the outflows and spin 
evolution of an accreting T Tauri star ([Pearson &: King 1995|) . 

Our calculation is, in principle, similar to the latter work because the gravitational effects 
of the companion star are removed from the code, however the purpose of this work and 
therefore the simulation details are quite different. We are interested in obtaining a steady- 
state structure of the disc itself and the spectral shape of the emitted radiation after the 
spin period of the central star has settled onto its equilibrium value. This calls for the 
determination of an appropriate spatial computational domain, a suitable particle injection 
scheme (initializing the calculation) and the need for a method to assess the temperature 
(and from it, the integrated spectrum) of the accreting matter. 

The inner and outer boundaries of our computational region are the stellar radius, i?*, 
and the ’’escape” radius, i?out (which is a parameter typically set at ~ 10i?co), respectively. 
Particles are injected into the computational domain (in a manner described below) and if 
during the course of the evolution they are found to leave this domain, they are removed 
from the calculation altogether. Once in the computational zone, each particle is subject 
to the held forces formulated explicitly in the rotating frame by the equations of motion 
^ through 0. In addition, short-range forces, rehecting the interaction between adjacent 
particles and thus simulating the hydrodynamical aspect of the dynamics, are employed 
using the HYDISC version of the particle hydrodynamics approach ( [Whitehurst 1988[ ). These 
forces are controlled by three parameters: C, determining the pressure force; Q, a constant 
of order unity and r^, the ’’chaining” mesh size. Physically C is the sound speed, Q is 
in fact the coefficient of restitution of the particle-particle collisions and is less than 1 for 
inelastic collisions, which must be assumed to simulate viscous flows, and is the viscous 
scale-length (for details, see Whitehurst 1988). 

The fact that HYDISC separates the treatment of the long-range (held) forces from 
the short-range (hydrodynamical) ones makes this code convenient for application to our 
problem. Most changes are introduced in the held part (as explained in the previous section 
on the equations of motion) and in the initialization schemes, that is, particle injection 
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(see below), while the treatment of particle interactions remains essentially the same as in 
HYDISC, with an appropriate choice of the parameters. 

The structure of HYDISC, as explained above, allows also to estimate, relatively easily, 
the energy dissipation (see below). 


3.2 Particle injection and disc initialization 


Our problem involves an accretion disc around a single star and thus the injection scheme 
has to be chosen accordingly. The disc around a T Tauri star is not being continuously 
replenished by an external stream, as is the case in discs fed by the secondary star in a close 
binary system, like a CV. This said, we still may chose between two distinct possibilities: 
a ’’burst” injection, in which all particles are injected simultaneously (or almost so) and 
a continuous injection. In both cases we have also to select the spatial distribution of the 
injected particles. Clearly, in order to obtain meaningful results the ultimate conhguration 
should not depend on the initial conditions. Computation time limitations dictate that the 
computational zone must be limited (as well as the total number of particles). Thus, the 
seemingly best option, i.e. to continuously inject particles at a radius much larger than i?co 
and wait for the disc to be built up by viscosity, is prohibitive. We have adopted instead 
what we call a ’’modihed burst” injection scheme. 

The particles were injected at an equatorial ring around a radius i?inj with a Gaussian 
spread in the individual particles injection radii and a small Gaussian displacement out of 
the disc plane as well. Each particle was initially assigned an azimuthal Keplerian velocity 
plus a random contribution, appropriate for the local disc temperature. The temperature 
of the disc, used to calculate the above mentioned random contributions to the injected 
particles’ velocities and the sound speed (and thus the constant C), was estimated using 
the analytical standard disc model. The actual formula which we have used is identical to 
formula 6 of PK: 


T(r) = T(3) 


27 


Ll-3-1/2 


r 3(1 - 


1/4 


(13) 


where r is the radius in units of i?* and T(3), is the temperature at 3 stellar radii, where 
true estimates are available ( |Basri fc Bertout 19891) . Although the temperature in our disc 
is somewhat different, this estimate actually influences only the pressure in the disc (in 
addition to the initial condition) and therefore the error it introduces should be negligible 
(see PK). 
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In order not to lose too quickly significant numbers of particles from the simulations, the 
’’modihed burst” scheme consisted of extending the burst duration to a hnite (typically a few 
stellar rotation periods) time Tb, with the injection rate decaying with time in a Gaussian 
manner. In addition, particles were added at a constant (but rather slow as compared to the 
burst) rate iVcont so as to compensate for particles lost from the computational zone (accreted 
or expelled). The spatial distribution of the added particles was kept hxed in time. Thus the 
modihed burst scheme actually gave rise to the following number of injected particles as a 
function of time (expressed here in units of the rotation period) 

iVinj(t) =fVb(t)+fVcontt, (14) 

where iVb(t) = /q iVo exp(—t^/2rb) dti . Typically, in our simulations the total number of 
particles in the burst (i.e. iVb(t = Sr?,) say) was between 5000 and 10000 and iVcont = 150 
particles per rotational period (see below). 

3.3 Calculation of the emitted spectrum 

With the present DBA scheme it is impossible to obtain a detailed thermal structure of the 
accretion how and therefore we aim at approximating just the ehect of the magnetic held on 
the emitted spectrum, relatively to the non-magnetic {k = 0) case. We assume that the disc 
is optically thick in the region of interest and that the energy dissipated in it by the viscous 
interaction is immediately locally radiated out in the vertical direction. These assumptions 
are essentially the same as in the standard Shakura & Sunyaev treatment of geometrically 
thin and optically thick accretion discs. What is diherent here is the fact that the viscous 
dissipation rate must be computed from the details of the numerical particle scheme and 
can not be found directly from the disc’s viscous how. As explained above, a part of the 
kinetic energy of the particles (blobs) is lost in collisions as these are inelastic (Q < !)• After 
dividing the computational zone into cells, by a two-dimensional grid in the disc plane, one 
can calculate the energy released in each cell AEij during a given time At, provided the 
masses of the particles are known. The composite spectrum can then be found by adding 
up the contributions of all the cells in the computational domain, with each cell assumed to 
contribute a black-body emission spectrum, corresponding to its ehective temperature 
„ _ ( A£)„ \ 

where Aij is the cell surface area and a is the Stefan radiation constant. 


( 15 ) 
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As explained in §2 (see equations and we have lumped the unknown values of the 
blob properties into a single parameter and chosen a specihc spatial dependence for it. Thus 
for a given value of the coefficient k, different combinations of the constituent parameters 
(and among them the blob’s mass) is allowed. Assuming some hxed average blob mass can 
thus give a spectral shape, which should be considered as a qualitative estimate and have 
a relative meaning, that is, the spectral shapes for different values of k can be perceived as 
an assessment of the effect of the magnetic held strength on the shape of the spectrum. In 
addition, since our computational zone does not include (due to computing time limitations) 
the extended outer portion of the disc, the results must be interpreted remembering this 
fact. In particular, since the outer cool parts of accretion discs contribute mainly to the long 
wavelength part of the emitted spectrum, our computed spectra clearly underestimate the 
spectrum strength for large A (see below in Fig. 3). 


4 THE SIMULATIONS 
4.1 Simulation parameters 

In choosing the parameters of the simulations we have tried to optimize the sometimes 
conhicting demands to obtain, on the one hand, meaningful results and on the other hand 
to limit, as much as possible, the CPU demands of the numerical calculations. 

We have set the stellar parameters in all the simulations to be M* = Mq, i?* = 1.5i?Q, 
F* = 10® s ~ 11 days, giving for the corotation radius Rco ~ 1.5 lO^^cm ~ 14F*. As explained 
before, the corotation radius serves as our length unit. 

Experimenting with different injection radii (the location of the centre of the injection 
ring), we have found that with the injection scheme (pT]) an injection radius of Finj = 1.4Fco 
is effective in building up rather quickly an accretion disc, which does not differ essentially 
from the one resulting (after a longer calculation) from a larger injection radius. In these 
tests we have shut off the magnetic interaction, but remembered from previous works that 
the injection radius should not be too close to corotation. Consequently in all the simulations 
we have used this value of the injection radius, with a Gaussian spread of up to iO.lFco- 

As discussed in the previous section, the parameters related to the particle scheme treat¬ 
ment of the pressure and viscosity include (7, Q and r^. We chose C to be equal to the 
isothermal sound speed, derived from the standard Shakura & Sunyaev disc temperature 
distribution. Q and r^ determine the viscous effects and have to be chosen with care. Q, the 
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restitution coefficient has a global effect on the simulation and since its reasonable values 
should lie between 0.5 and 1 (Q < 0.5 allows particle interpenetration and Q = 1 corre¬ 
sponds to elastic collisions). We have found that changing Q within the above mentioned 
range has only a small effect on the final quasi steady state, but it influences the length of 
the time until such a state is reached. Whitehurst (1988) reached a similar conclusion by 
experimenting with Q = 0.6 and 0.7 values in the original HYDISC code. Thus, wishing to 
shorten the simulation time, without, however, making the viscosity unphysically large, we 
have opted for Q = 0.7 and used this value in all our simulations, save for a few test cases. 

The parameter r^, setting the viscous interaction length scale, has a local effect on 
the simulation. When deciding on the range of variation for this parameter we have taken 
into account the fact that too small a value of causes the viscous evolution time to be 
unacceptably long. In addition, because the number of cells in the computational range 
becomes large, one has to increase the number of particles in order to keep the particle 
treatment meaningful. Both effects increase the CPU time and we have found that for 
Tv ^ 0.02 we were unable to efficiently compute the disc evolution. On the other hand for 
approaching the value of 0.1 the increase in viscosity causes the particles to clump together 
and the evolution becomes unrealistically viscosity dominated. Performing test calculations 
with no magnetic interaction {k = 0) we found that an initial ring injected around 1.4i?co 
viscously spreads to a structure very similar to a standard thin disc in a time of ~ 100 — 300 
rotational periods when is kept within the range 0.08 >r^> 0.03. Using such a calculation 
it is also possible to find what is the effective a (the viscosity coefficient in the Shakura & 
Sunyaev prescription) with a given choice of Q and r^. A test calculation with = 0.05 (in 
units of Rco, see below) and Q = 0.7 was very similar to the evolution of a ring of matter to 
as it is found from a standard exercise with an a viscosity (see e.g. Pringle 1981) value of 
~ 0.1. 

The parameters k and were varied in the different simulations. The range of the 
variation of k was chosen so as to correspond to polar stellar magnetic value of up to several 
kiloGauss (when the other parameters in k are set to their representative values, see equation 
1^ and Tv was chosen within the range mentioned above. 
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4.2 Results of representative cases 

We have performed 12 full simulations of an accretion disc in an aligned dipole magnetic field 
for different choices of the pair of parameters k and r^. The simulations were initialized using 
the injection scheme described in equation |T^ with parameters chosen so as to insure that the 
number of particles does not drop below ~ lO'^ during the simulation. The simulations were 
continued until a quasi-steady configuration was achieved, that is, no significant changes of 
the structure were observed over a significant time (see below). We have chosen to describe 
here in some detail the simulations with = 0.05 (simulation group 3 in our list), that is, 
with the viscosity strength chosen in the middle of its feasible range. As discussed in the 
previous section this choice corresponds to an effective a ~ 0.1. Three runs were performed 
for different values of the magnetic interaction strength [k = 0,10,100), where the non¬ 
magnetic case {k = 0) was run to serve as reference, in particular for the spectra, which 
have quantitative meaning only in comparison to each other (as explained before). 

These three runs were all allowed to continue for a total simulation time of 250P*. A 
quasi-steady state was achieved already after approximately half of this time. The CPU 
time requirement for each run, with the number of particles kept at ~ 10"^ (the accreted and 
ejected particles were replenished as explained above), was about 100 hours on a Sun ULTRA 
2 workstation, when the program was coded in C, to increase its efficiency. In Fig. 1 we show 
the surface number density of the particles (a quantity proportional to the disc surface mass 
density, customarily called S) as a function of radius for the three values of k. Note that in 
the innermost and outermost parts of the disc, regions in which the surface density curves 
are dashed, the number of particles is too small to render the result to be quantitatively 
reliable. We also plot on this graph the analytical result, arising from a standard thin disc 
formula with inner cutoff (see below). This is represented by the fully dashed curve, which 
which follows the k = 100 case for most of the range. 

From this figure it is evident that the effect of the magnetic field on the disc structure 
is important, as it significantly modifies the functional dependence of the surface density 
on the radius. For fc = 0 (no magnetic field) the particles form a rather ’’collapsed” with 
its density growing sharply towards the central object. This is so because our simulations 
were tailored for the high k case and therefore the particle injection rate was too small 
for the buildup of a steady disc in a small k case. A growing magnetic field strength has a 
dynamical effect on the disc. As could be expected from a previous study, in which individual 


14 Ultchin et al. 



Figure 1. Surface number density as a function of radius (in units of -Rco) for rv = 0.05 (simulation 3) and three values of 
k. The number of particles in the innermost and outermost regions, denoted by dashed portions of the curves, is too small to 
render the results to be quantitatively reliable there. The fully dashed curve is an analytical result of a disc truncated at Rco 


orbits of particles were calculated ([Ultchin et al. 1997|) , and from the DBA simulations of 
PK, the magnetic interaction drags blobs residing inside corotation towards the star and 
pushes out blobs whose orbits are outside corotation. These effects are modified here by the 
fact that the dipole magnetic field decays outwards and by the interblob interactions. The 
combined effect is a significant depletion of blobs close to the star (a ” hole”) and an increase 
of density at large radii (a ’’flattening” of the functional dependence of the surface density 
on radius). The size and depth of the hole and the extent of the flattening of the surface 
density obviously depend on fc, That is, on a combination of the disc density (which depends 
on the accretion rate) and the strength of the magnetic field. 

Analytic estimates of the hole radii at spin equilibrium, as given by Cameron and Camp¬ 
bell (1993), predict that the disc should end abruptly close to the corotation radius (but 
depending also on the fastness parameter of the star). Our numerical results do not yield an 
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R/Rco 


Figure 2. The effective temperature (in arbitrary units) as a function of radius (in units of -Rco) for = 0.05 (simulation 3) 
and k = 100. The dashed curve is the analytically obtained temperature distribution of a classical accretion disc, cut off inside 
at = Rco (see 0 ) 


abrupt disc termination, but the value of ~ Rco is certainly consistent with these results. In 
the following comparisons of our results to analytical solutions we have used this value of 
the inner cutoff radius, that is, the analytical formulae giving rise to the fully dashed curves 
in Figs. 1 and 2 (see below) result from 


m) 



( 16 ) 


and 


T^{R) = 


3GMM 


\J Rco/R 


(17) 


‘SnRfa 

with the values of v and M chosen appropriately, from our calculation, so as to enable the 
comparison. 

If Fig. 2 we show the comparison of the analytical (effective) temperature distribution 
as given in (p!7|) with the numerical result for the case of fc = 100, when the outer parts 
of the disc resemble most the analytical result (see Fig. 1). The lack of the inner, hottest. 
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part in the numerical results is apparent. In these calculations the method described in §3.3 
was used to monitor the energy dissipation during the simulations and calculate the discs 
temperature structures and their spectra. The spectral flux distributions {\F\) in the three 
cases of k depicted in Fig. 1 are displayed in Fig. 3. 

It is apparent from the figure that with increasing k the spectral distribution of the 
emitted radiation is shifted towards longer wavelengths. This trend is quite clear since as 
we have seen in Fig. 1 and 2, the hottest material, which resides close to the star, is sig¬ 
nificantly depleted by the magnetic interaction. In addition, as a result of the increase of 
the surface density in the outer cooler parts, the contribution to the long-wavelength part 
of the spectrum is more prominent. We stress once again that the temperature calculation 
and hence that of the spectrum suffers from the fact that the mass of the blobs is essentially 
unknown, so that the results displayed in Fig. 3 are of a qualitative nature and their value 
is mainly in uncovering the trend of the spectral behaviour. In addition, computer time 
limitations do not allow us to extend the calculation to the outer parts of the disc (here, in 
simulation 3 we have i?out = 10i?co)- Thus the spectral energy distributions in Fig. 3 fail to 
account for the cool material lying outside the outer limit of the computational zone and 
which should have an effect of increasing the strength at long wavelengths. Since the disc is 
more extended for large k this effect should matter the most in the strongly magnetic case. 
The actual structure of the accretion disc for the three values of k can be seen in Fig. 4, in 
which maps of the the surface density and temperature are shown. 

From the examination of Fig. 3 and the above considerations it is thus reasonable to con¬ 
clude that the magnetic interaction, at least as it is dealt with in this work, tends to flatten 
the spectrum for long wavelengths, consistently with the observed IR excess mentioned in 
the Introduction. 

The results of all the other simulations (i.e. the ones with a different value of the viscosity 
scale parameter are similar to the ones described above and thus the conclusions are 
similar. The actual value of has only a little effect on the final (i.e. quasi-steady) surface 
density and temperature distribution, as long as this value is within the calculation feasibility 
range (see above). The difference is in the time that has to pass until the configuration settles 
down. This is understandable as the evolution time is essentially the viscous time scale and 
the latter decreases with increasing r^. 
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Figure 3. Spectra of the accretion discs in simulation 3. AF(A) in arbitrary units (see text) is plotted as a function of wavelength 
A in microns for three values of the magnetic interaction parameter k. 


4.3 Effects of special magnetic field geometry 

We have performed, in addition to the simulations described above in which the magnetic 
held was that of an aligned dipole, two similations with different magnetic conhgurations. 
In the hrst one we have used a magnetic dipole whose axis was tilted to the disc axis and in 
the second one a ’’dipole slice” (see below) was introduced in order to estimate the effects 
of a localised magnetic held conhguration. 

4.3.1 Tilted dipole 

A single simulation was performed with = 0.05; k = 100, that is, similar to one of the 
simulations described above, but with the magnetic dipole axis tilted at an angle of 45° to 
the disc axis. The results show relatively few changes as compared to the aligned dipole case. 
There was almost no perturbation of the disc matter, relatively to the aligned dipole case, 
far out of the corotations radius. This is understandable since the blobs orbiting at these 
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Figure 4. Surface density (left panels) and temperature (right panels) maps for simulation 3 (rv = 0.05) for the three values of 
the magnetic interaction strength. Arbitrary units are used and in the surface density maps a side view of the disc is also given, 
in the upper part of each panel.Surface density (left panels) and temperature (right panels) maps for simulation 3 (rv = 0.05) 
for the three values of the magnetic interaction strength. Arbitrary units are used and in the surface density maps a side view 
of the disc is also given, in the upper part of each panel. 
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positions experience in this case a magnetic field whose direction oscillates with the period of 
the stellar rotation (which is considerably shorter than the blob’s Keplerian period). Tims, 
it is expected that the blobs actnally react to the average magnetic field, i.e. behave in a 
similar way as in the aligned dipole case. 

This expectation is wrong for blobs close to and inside corotation. Indeed, their motion 
shonld be channeled there by the magnetic field. This effect conld not be observed in onr 
simnlations since the density of the particles in the magnetosphere, above the disc plane, is 
very small. A very high nnmber of particles or alternatively an extremely strong magnetic 
field is needed in order to resolve, in onr type of calcnlation, the particle flow along field 
lines. This remark obvionsly applies to the aligned dipole case as well. 

4-3.2 Dipole slice 

To simnlate the case in which the actnal stellar magnetic field resembles a localized loop, 
similar to the case stndied by PK, we have nsed a ’’slice” of the aligned dipole, with Ganssian 
fall-off in strength in both the radial coordinate R and the azimnthal angle 6 . Thns we have 
performed a simulation with = 0.05; k = 500 and with the magnetic field centered on the 
position Rq, 9q 

Bsiice = 5dipoieexp[-(6' - 9of/2dl] exp[-(i? - Rof/2dl], (18) 

where Sdipoie is a dipole field (see the functional dependence in equation ^ and dn, de are 
the radial and angular Gaussian fall-off scales in the radial and azimuthal directions. In the 
simulation we have chosen i?o ~ Rco, consistently with the findings of KR that this should 
be the equilibrium configuration. 

The results clearly depend on dji and d^. When we chose ~ Rco, the disc beyond Rco 
resembled the non-magnetic case (A; = 0), since the field is practically zero there. Gloser to 
corotation and inside it, the density distribution looked like that of a full but weaker dipole, 
since the dipole slice field effectively acts on the orbiting blobs only along a fraction ~ de/2TT 
of the blob’s orbit. 

5 SUMMARY AND CONCLUSIONS 

The most prominent feature of our simulation results is the formation of a low density region 
(a hole) in the inner part of the disc as a result of the magnetic interaction. In addition, 
the density distribution becomes more spread out, as the slope of logS(R) decreases with 
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increasing k. This behaviour is typical and essentially qualitatively independent of all the 
other parameters. Although the ’’hole” is the more noticeable feature (see e.g. the density 
maps in Fig. 4), it is less signihcant to overall observable features than the global change in 
the density distribution. 

We have also seen, by comparing our calculations with what can be expected from 
analytical results, that we uncover some features (like the absence of hot material close to 
the inner hole) which can not be obtained from just classical disc models with a hole inside. 

Due to computer power restrictions we were unable to achieve a high enough resolution, 
so as to see the accretion flow along the held lines and outhows from the system. In addition 
we were able to simulate only a limited portion of the disc and thus the calculations of 
the spectrum reveal only the general trend. This trend, resulting from the battening of the 
surface density distribution, was always to shift emission power toward longer wavelengths. 
The battening of XFx is most apparent near its maximum and it is reasonable that were 
it not for the close (computational) cutob in R, the shape would remain bat to longer 
wavelengths. We propose therefore that the IR excess in T Tauri stars can be attributed to 
magnetic interaction, which modibes the functional dependence of the surface density in the 
surrounding discs. 

Existing models attempting to explain the IR excess in T Tauri stars fall into two dis¬ 
tinct categories: those invoking geometrical factors and others, proposing energy dissipation 
mechanisms operating preferably in the outer parts of the disc. Our model suggests a cor¬ 
relation between the spectrum batness and the magnetic held strength, appears to be quite 
robust (practically any shape of the magnetic held would do) and does not require assump¬ 
tions about bared shapes of discs or unusual energy dissipation modes. At this stage, the 
results of our calculation provide little more than support for a qualitative promising idea. 
As it was mentioned in the Introduction, it is possible to apply a similar prescription to 
the case in which the magnetic held penetrates the disc. We can reasonably expect that the 
results of such a calculation should not be too diberent, at least qualitatively, from the ones 
presented in this paper. It appears that all that is required is magnetic held lines imparting 
a torque on the gas, which changes in sign as we cross the corotation radius. 

Signihcant results for both models of the magnetic interaction can be achieved in high res¬ 
olution (signibcantly larger number of particles). However, to extend the idea into a reliable 
quantitative model, full multidimensional MHD simulations, including radiative transfer, 
have to be ultimately performed. 
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